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Abstract 

This work continues a program to systematically generalize the Skyrme Hartree-Fock method for 
medium and heavy nuclei by applying effective field theory (EFT) methods to Kohn-Sham density 
functional theory (DFT). When conventional Kohn-Sham DFT for Coulomb systems is extended 
beyond the local density approximation, the kinetic energy density r is sometimes included in 
energy functionals in addition to the fermion density. However, a local (semi-classical) expansion 
of T is used to write the energy as a functional of the density alone, in contrast to the Skyrme 
approach. The difference is manifested in different single-particle equations, which in the Skyrme 
case include a spatially varying effective mass. Here we show how to generalize the EFT framework 
for DFT derived previously to reconcile these approaches. A dilute gas of fermions with short-range 
interactions confined by an external potential serves as a model system for comparisons and for 
testing power-counting estimates of new contributions to the energy functional. 
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I. INTRODUCTION 



Density functional theory (DFT) is widely used in many-body applications with Coulomb 
interactions because it provides a useful balance between accuracy and computational cost, 
allowing large systems to be treated in a simple self-consistent manner. In DFT, the particle 
density p(x) plays a central role rather than the many-body wave function. DFT has the 
generality to deal with any interaction but has had little explicit impact on nuclear structure 
phenomenology so far (see, however Refs. [l[ 0, SB), although the Skyrme-Hartree-Fock 
formalism is often considered to be a form of DFT [5|. In previous work, effective field theory 
(EFT) methods were applied in a DFT framework, with the ultimate goal of calculating 
bulk observables for medium to heavy nuclei in a systematic fashion. The present work 
takes another step toward this goal by generalizing the EFT framework for DFT to include 
the kinetic energy density in the same way it appears in the Skyrme approach. 

Kohn-Sham (KS) DFT maps an interacting many-body system to a much easier-to-solve 
non-interacting system. Hohenberg and Kohn proved that the ground-state density p(x) 
of a bound system of interacting particles in some external potential f (x) determines this 
potential uniquely up to an additive constant 0|. They showed that there exists an energy 
functional E[p\ that can be decomposed as 

E[p{^)] = Fhk[p(x)] + j rf^xt;(x)p(x) , (1) 

where the functional -Fhk[p] is called the HK free energy, and is universal in the sense that 
it has no explicit dependence on the potential f (x). A variational principle, formulated 
in terms of trial densities rather than trial wavef unctions, ensures that the functional E[p\ 
is a minimum equal to the ground-state energy when evaluated at the exact ground-state 
density. 

The free-energy functional -Fhk[p] can be further decomposed into a non-interacting ki- 
netic energy Ts[p\, the Hartree term Eh[p\, and the exchange-correlation energy Exc[p\ 0]: 

Fhk [p] = [p] + Eh [p] + E^, [p] . (2) 

Kohn-Sham DFT 0, Hll allows Ts[p\ to be calculated exactly [see Eq. (jHUj) ]. which 

efficiently treats a leading source of non-locality in the energy functional. Given the Hartree 
and exchange-correlation functionals, the solution of the interacting problem at zero tem- 
perature reduces to solving a single-particle Schrodinger equation, 

+ Vsij^)] ^/3(x) = ep'ippipi) , (3) 



2M 



for the lowest A orbitals ^/^(x) (including degeneracies),^ where the effective local external 
potential fs(x) is given by 

Vs{^) - nxj + ^^^^^ + ^^^^^ . i4j 



^ This is only true in the absence of pairing. See Ref. for a discussion on generahzing DFT/EFT to 
accommodate pairing correlations. 
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The wave functions for the occupied states in this fictitious external potential fs(x) generate 
the same physical density p(x) through p(x) = l'^/9(^)P ^^^^ original, fully 

interacting system in the external potential f (x). 

The practical problem of DFT is finding a useful explicit expression forE,,[p] BiHini. 

For Coulomb systems, the conventional procedure is to approximate E^c in the local density 
approximation (LDA) by taking the exchange-correlation energy density at each point in 
the system to be equal to its value for a uniform interacting system at the local density 
(which is calculated numerically) 01 , and then to include semi-phenomenological gradient 
corrections 0, These gradient corrections have become steadily more sophisticated. 
The initial attempts at gradient corrections violated a sum rule p, lul , which was fixed by 
the generalized gradient expansion approximation (GGA), where spatial variations of p(x) 
are constrained by construction to conform with the sum rule To further improve the 
exchange-correlation functional, the semi-local kinetic energy density t(x) = |V'?/'/3(x)p 
was built into the GGA formalism to construct meta-generalized gradient approximations 
(Meta-GGA) [l^. However, in practice, the kinetic ener gy d ensity is replaced by its local 
expansion in terms of the density p(x) and its gradients so that the energy is treated 
as a functional of the density alone. The Kohn-Sham single-particle equation still takes the 
form of Eq. Q . 

The Skyrme-Hartree-Fock approach to nuclei is also based on energy functionals and 
single-particle equations, which suggests a link between the traditional DFT and nuclear 
mean- field approaches p, U^, [l^ • In calculations with the usual density-dependent Skyrme 
force, the energy density for spherical, even-even N = Z nuclei takes the form of a local 
expansion in density |l8|, |l9| , 



— r(x) + -tolpU)? + — t3[p(x)]'+" 
2M ^ ' 8 ^ ^ 16 ^ ^ 



— (3ti + 5t2)p(x)r(x) 



16 



-l(9ti - 5t2)|Vp(x)p - ^W-op(x)V ■ J(x) + - t2)[J(x)] 



(5) 



(see Ref. ji^l for a more general treatment). The density p, kinetic density r, and the 
spin-orbit density J are expressed as sums over single-particle orbitals '?/'/3(x): 



p(x) ^ i^/^w 



r x 



/3 



J(x)^5^^t(^)(Vxa)^^(x), (6) 



where the sums are over occupied states. The tj's, Wq, and a are generally obtained from 
numerical fits to experimental data. Varying the energy with respect to the wavefunctions 
leads to a Schrodinger-type equation with a position-dependent mass term [iM Eoj]: 



1 



2M*fx 



-V + f/(x) + l^(x)U^(x) = £^z^^(x 



(7) 



where M*(x) is given by 



2M* 



2M 



+ 



3 5 

— ti H t2 

16 16 



P(x) 



(8) 



and W{x.) is a spin-orbit potential (see Ref. j20| for details). The appearance of M*(x) and 
W{x) are a consequence of not expanding r and J in terms of p. 
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The Skyrme approach has had many phenomenological successes over the last 30 years 
and, generahzed to include the effects of pairing correlations, continues to be a major tool 
for analyzing the nuclear structure of medium and heavy nuclei [21„ ,2^ ,2^ 2^ ■ The form of 



the Skyrme interaction was originally motivated as an expansion of an effective interaction 
(G-matrix) in the medium [2^. Negele and Vautherin made the connection concrete with 
the density matrix expansion method [2^, but there has been little further development 
since their work. Many unresolved questions remain, which become acute as one attempts 
to reliably extrapolate away from well-calibrated nuclei and to connect to modern treatments 
of the few-nucleon problem. Why are certain terms included in the energy functional and 
not others? What is the expansion parameter (s)? How can we estimate uncertainties in the 
predictions? Are long-range effects included adequately? Can we really include correlations 
beyond mean-field, as implied by the DFT formalism? We turn to effective field theory to 
address these questions, using a DFT rather than a G-matrix approach as in Ref. (27| . 

Effective field theory (EFT) promises a model-independent framework for analyzing low- 
energy phenomena with reliable error estimates (isl . I29I 30l . 3ll Is^ . In Ref. [s^ , an effective 



action framework was used to merge effective field theory with DFT for a dilute gas of 
identical fermions confined in an external potential with short-range, spin-independent in- 
teractions (see also Refs. [s^,!!^). The calculations in Ref. 0| used the results from the 
EFT treatment of a uniform system in Ref. 3^| in the local density approximation. In the 
present work, we extend the formalism to include the kinetic energy density as a functional 
variable, leading to Kohn-Sham equations with position-dependent M*'s (see Ref. ||37| for 
an earlier discussion of such an extension to Kohn-Sham DFT). We take a (small) step 
beyond the LDA by evaluating the full r dependence in the Hartree-Fock diagrams with 
two-derivative vertices, which leads to an energy expression similar to the standard Skyrme 
energy (excepting the spin-orbit parts). 

The plan of the paper is as follows. In Sect. ^1 we extend the EFT/DFT construction 
for a dilute Fermi system to include a source coupled to the kinetic energy density operator. 
A double Legendre transformation, carried out via the inversion method, yields an energy 
functional of p and r, and Kohn-Sham equations with M*(x). In Sect. IIIIl we illustrate the 
formalism with numerical calculations for a dilute Fermi system in a trap, including power- 
counting estimates. Section HVI is a summary. For the bulk of the paper, we will restrict our 
discussion to spin-independent interactions but in Sec. IllDI the extension to more general 
forces is discussed. 



II. EFT/DFT WITH THE KINETIC ENERGY DENSITY 

A dilute Fermi system with short-range interactions is an ideal test laboratory for effec- 
tive field theory at finite density, but it is also directly relevant for comparison to Skyrme 
functionals. We describe such a system using a general local Lagrangian for a nonrela- 
tivistic fermion field with spin- independent, short-range interactions that is invariant under 
Galilean, parity, and time- reversal transformations j32L 1361]: 



2MJ 



Co 



G 



16 



^-Tr(^V)^ + TF (^^)"^(^VV)+ h.c. 



(9) 
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where V is the Gahlean invariant derivative and h.c. denotes the Hermitian con- 

jugate. The terms proportional to C2 and C2 contribute to s-wave and p-wave scattering, 
respectively, while the dots represent terms with more derivatives and/or more fields. To 
describe trapped fermions, we add to the Lagrangian a term for an external confining po- 
tential f (x) coupled to the density operator v^^ip^ifj js^. For the numerical calculations, 
we take the potential to be an isotropic harmonic confining potential, 

t;(x) = ^mtu^ |x|2 , (10) 

although the discussion holds for a general non-vanishing external potential. This La- 
grangian can be simply generalized to include spin-dependent interactions, as described in 
Sec. Ill D| although the complete set of terms grows rapidly with the number of derivatives. 



A. Effective Action and the Inversion Method 

We introduce a generating functional in the path integral formulation with C from Eqs. Q 
and (fTUj) supplemented by c-number sources J(x) and ?7(x), coupled to the composite density 
operator and to the kinetic energy density operator, respectively. 

For simplicity, normalization factors are considered to be implicit in the functional integra- 
tion measure 3^ 3^. The fermion density in the presence of the sources is 

p[x)^{i;\x)i;{x))j^,= ^-^^^ . (12) 
and the kinetic energy density is given by 

r{x) ^ (V^t(^) . VV^(x)),,, = ^^1^ , (13) 

where functional derivatives are taken keeping the other source fixed. The effective action 
is defined through the functional Legendre transformation 

r|p..]^Hai.,]-/A.(.M.)-/A»,. (14) 

which implies that F has no explicit dependence on J and 77. 

As in Ref. ji^, we choose finite-density boundary conditions that enforce a given particle 
number A by hand, so that F and W are functions of A and variations over p{x) and t{x) 
conserve A. This is naturally achieved by working with a fixed number of Kohn-Sham 
orbitals, as seen below. By limiting ourselves to time independent sources and densities, we 
can factor out a ubiquitous time factor corresponding to the time interval over which the 
source is acting, and write ji^ 



F[p,r] = F[p,r]x 



00 

dt 

■00 



-1 



-E[p,T]. (15) 
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and similarly with 1^[J, r^] and the expansions below. The effective action has extrema at 
the possible quantum ground states of the system, and when evaluated at the minimum is 
proportional (at zero temperature) to the ground-state energy jsl, |3^ . In particular, 
Eq. (fl3|) defines an energy functional E[p, r], which when evaluated with the exact ground- 
state density p and kinetic energy density r, is equal to the ground-state energy. 

If we take functional derivatives of Eq. (jl4p with respect to J(x) and 77(x), and then 
apply Eqs. (|12p and (jl3|) [s^, we obtain two equations that can be written in matrix form 
as 



/ 5p{y) 6r{y) \ 



5J(x) 
\ 5r/(x) 



(5J(x) 
5r/(x) / 



5p(y) 



^(y) 



6T\p,t] I . 



(16) 



The invertibility of the transformation from {J, to {p, t} implies that the matrix in 
Eq. (fT^ has no zero eigenvalues, which means that the elements of the vector vanish iden- 
tically for each x, so that 



5p(x) 



-J(x) 



(5rfx) 



-77(x) . 



(17) 



Thus, the effective action when evaluated at the exact ground state density and kinetic 
energy density is an extremum when the sources are set to zero, which corresponds to the 
original (source-free) system. The convexity of F implies that the energy (equal to minus 
r[p, r] at the extremum) is a minimum.^ 

These properties are analogous to conventional applications of effective actions where the 
Leg endre transformation is with respect to one of the fields rather than a composite operator 
0,0,0]. A possible complication in the present case would be if new divergences arose 
in ly [J, 77] , which generally happens when adding a source coupled to a composite operator 
lisi 146^ . However, that is not the case here [12]. On the other hand, sources coupled to 
operators such as tpil) -|- tp'^il)'^ ^ which arise when considering pairing, will introduce new 
divergences, including for the kinetic energy density. This issue will be considered elsewhere 



To carry out the inversion and Leg endre transformation, we apply the inversion method 



of Fukuda et al. [38|, |39|, |41|, |47|5 |48|, |49[. This method applies to any system that can 



be characterized by a hierarchy, which we organize by introducing a parameter A that is 
ultimately set to unity. This parameter can label different orders in a coupling constant 
expansion (e.g., powers of for the Coulomb interaction), in a large iV expansion, or in an 
effective field theory expansion. Here we apply the latter, and associate powers of A with 
the orders in the dilute EFT expansion derived in Refs. [s^ and [s^. The effective action 
is given a dependence on A 

f = f[p,r,A], (18) 
which is treated as an independent variable. The Legendre transformation defining F follows 



^ In Ref . |41| , a proof is given of the invertibility of the Legendre transformation for the EucHdean version 
of the functions. The proof extends directly to our Minkowski functions with any number of sources, as 
long as they are coupled linearly. 
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LO: {^^^"C^ NLO: 
(a) 

NNLO: 0000+ 




+ 




FIG. 1: Hugenholtz diagrams for a dilute Fermi gas through next-to-next-to-leading-order (NNLO) 
in an EFT expansion. 



from Eq. dH)): 



r[p, r, A] = W[J, 77, A] - / rf^x J(x) p(x) - / rf^x r/(x) r(x 



(19) 



where J and 77 depend on A as well as being functionals of p and r. 

Now we expand each of the quantities that depend on A in Eq. (fT^ in a series in A |3^ 
treating p and r as order unity, and substitute the expansion for J and rj into the expansion 
for W. Equating terms with equal powers of A on both sides of Eq. (|19p after carrying out 
a functional Taylor expansion of W[J, rj] about Jq and rjo gives a series of equations relating 
the Ti, Wi, Ji, and r]i, where the subscript / indicates the power of A. These equations allow 
the r^'s to be constructed recursively (see Ref. for explicit expressions). Since p and r 
are independent of A, the sources for any I satisfy 



5ri[p,r] 
" 5p(x) 



STi[p, 



T 



(5r(x) 



(20) 



which is just the term-by-term expansion of Eq. ()17|). All of the J^'s and 77; 's as defined 
here are functionals of p and r, and this functional dependence will be understood even if 
not explicitly shown from now on (and functional derivatives with respect to p will have 
r held fixed, and vice versa). We identify the Wi^s for Z > 1 in the present case with the 
diagrammatic expansion in Fig. [T] [s^, S]- That is, Wi is given by diagram (a), W2 by 
diagrams (b) and (c), and W3 by diagrams (d) through (j). 
The zeroth-order equation from Eq. ()19|) is 



To[p, t] = Wo[Jo, r]o] - / d^^ Jo(x) p(x) - / rf^x r]o{x.) r(x 



(21) 
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The corresponding zeroth order expansion of Eqs. ()12|) and ()13p is 



P(x) 



r X 



%(x) 



(22) 



(this also follows by taking functional derivatives of Eq. (PT|) with respect to p and r and using 
invertibility again, see [3Jl|). Note that these are the exact ground-state densities. Because p 
and r are treated as order unity, this is the only equation in the expansion of Eqs. ()12|) and 
(fT!?jl in which they appear. Thus, the sources Jo(x) and ?7o(x) are particular functions that 
generate the expectation values p and r from the noninteracting system defined by A = 0. 
(The existence of Jo(x) and rjolx.) is the heart of the Kohn-Sham approach.) The inversion 
method achieves this end by construction. 

The exponent in the non- interacting generating functional Zo[J, r]] is quadratic in the 
fermion fields. 



2M 



t;(x) + Jo(x) 



'^'^2M~ ^(x) + ^o(x) - V • r7o(x)V 



i!+ i d'^X T]Q{x.)Vi!^ ■ Vip 



(23) 



which leads us to define the Green's function G^^ of the Kohn-Sham non-interacting system. 
This Green's function satisfies 



idt + V ■ 



1 



2M*('x) 



-V - vi^ 



+ Jo(x)^ ^^.(xt, xY) = 53(x - x')5(t - t') 



(24) 



with finite density boundary conditions j50[ and a position-dependent effective mass defined 
by 

^ ^ -^o(x). (25) 



2M*(x; 

Kohn-Sham orbitals arise as solutions to 

1 



2M 



hV 



V + t;(x) - Jo(x)] ^/'fc(x) = Sk^ki?^) 



(26) 



where the index k represents all quantum numbers except for the spin. Note that Eq. (|26|) 
is in the form of the Skyrme single-particle equation [Eq. ((7j)] (without the spin-orbit part). 
The spectral decomposition of in terms of the Kohn-Sham orbitals is [3^ 



zG°,(xt, xY) = M^) ^ki^') e-*^*(*-*') m - t') Oisk - £f) - eit' - t) 9isF - e,)] • (27) 



It follows that VFo[<^0) Vo] Tr ln(G'°g)~^ (since it is quadratic it yields a simple determinant), 
corresponding to the system without interactions, and can be written explicitly in terms of 
the single-particle Kohn-Sham eigenvalues as 



Wo[Jo,Vo\ 



(28) 
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where v is the spin-isospin degeneracy, as expected for a system without interactions. Equa- 
tion ()22j) apphed to Eq. ()28|) [with the help of Eq. ()26|)] imphes that p(x) and r(x) follow 
as in Eq. © from the orbitals |^]. By using Eq. ()28|) in Eq. ()2H) and then eliminating 
using Eq. the lowest order effective action can be written two ways, 

occ. „ » 

fo[p,r] = -v^Y^Sk- / rf^x Jo(x) p(x) - / rf^x r7o(x) r(x) 

= -Ts\t\- y^^x t;(x) p(x) , (29) 



where 



^x rfx) 



is the total kinetic energy of the KS non-interacting system. 
The first-order equation in A from Eq. (jl9j) is 



(30) 



- y d^x Ji(x) p(x) - y (i^x ?7i(x) r(x) 
= W^i[Jo,r7o] . (31) 



The complete cancellation of the J\ and ri\ terms from applying Eq. ()22|) occurs for the J/ 
and ?7/ terms in the equation for Yi for all / > 1. Thus for a given /, we only need VT^'s with 
k less than or equal to /, and J^'s and r/^'s with /c smaller than /. The VF^ functionals are 
constructed using conventional Feynman rules in position space with factors as in Ref. js^, 
but with fermion lines representing iG\^. 

In particular, the LO effective action is given by [s^ 

fi[p,r] = l^i[Jo,r7o] = ^^^(^^-l)Co ^ t^^x G°,(x, x+) ^".(x, x+) . (32) 

The density can be directly expressed in terms of the Kohn-Sham Green's function with 
equal arguments as 

p(x) = -ivGlJyX.x"^) , (33) 

so that we have 

fi[p,r] = -^^^ Co 1 d^^ |p(x)p = flip] , (34) 
which is minus the Hartree-Fock energy. Using Eq. (j^Uj) . we obtain 

(x) = ^M^izll p(x) , (x) = . (35) 
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After canceling the J2 and 772 terms as advertised above, the second-order effective action 
is given by 



m[Jo,vo] 



+ 1 I d'^d'y ^\ 7_";7; ''^[^ Ji(x) Ji(y) + ; U'^d'y 



5Jo(x) 5Jo(y) 



^i(x) r/i(y) 



+ / d^x dy 



(5Jo(x) %(y) 

W^2[<^0; ''?o] is calculated from the graphs Figs. [Tfb) and (c): 

W2[Jo,Vo] = Jd'x d'y Gl{x,y)Gl{x,y)Gl{y,x)Gl{y,x) 

C2 



(36) 



- ti^ii, - 1)' ^ y A G°,(x,x+)G°,(x,i/X(y,x)GL(y,y+) . (37) 
The other terms in r2 completely cancel against the "anomalous" graph of Fig. ^c) so that^ 
r2[p,r] = - 1) ^ 1 d'x dS Gl{x,y)Gl{x,y)GUy,x)Gl{y,x) , (38) 

which is equal to the contribution made by the "beachball" diagram [Fig. CJb)] to the energy 
[up to the time factor of Eq. This cancellation is proven exactly as in Ref. [33] after 

generalizing to the matrix notation of Eq. Calculation of the third-order effective 

action in the inversion method similarly leads to cancellation of the "anomalous" graphs in 
W3 given by Figs. HJ^d), (e), and (f), leaving only Fig. mg) through (j) as contributors. All 
higher orders in r[p, r, A] are determined in a similar manner, as described in Refs. 

Bull- 

In order to solve for the orbitals in Eq. (pUj) and to calculate the energy, we need expres- 
sions for Jo(x) and 77o(x). Since J(x) = rj{y.) = in the ground state, Eq. (fTTj) becomes 
a variational principle that, together with Eq. (pUj). yields self-consistent expressions for Jq 
and 770 [33^ : 



Jo(x) 



^oixj 



1>1 



(^Fint [p,r] 



gs (5p(x) 
(^Fint [p,r] 



5r(x) 



(39) 
(40) 



gs 



where Fint[p, t] is the interaction effective action 

fint[p,^] = ^fi[p,T] 



(41) 



i>i 



^ This complete cancellation does not occur for long-range forces, or if the zero-range delta functions at the 
Co vertices are regulated by a cutoff rather than by dimensional regularization, as used here. 
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and the subscript "gs" refers to the ground-state. (Note that at this stage A = 1, and 
Eqs. fj39|) and ()40j) mix all orders of the original inversion-method expansion into the Kohn- 
Sham potentials Jo and ?7o.) A given approximation corresponds to truncating Eq. (PT|) at 
^max and then carrying out the self-consistent calculation. We refer to Zmax = 1 as leading 
order, or LO, /max = 2 as next-to-leading order or NLO, and /max = 3 as NNLO. 



B. EFT for Dilute Fermi Systems 

The non-interacting energy density at zero temperature for A particles with spin- 
degeneracy u in volume V can be written as 

where the density p is 



5 V2M 



For a uniform system, the order-by-order in kptts corrections to Eq. (j42j) due to interactions 
can be calculated in the EFT via the Hugenholtz diagrams in Fig.Hl as described in Ref. |36l | . 

For a dilute Fermi system, the coefficients Co, C2, and C2 can be expressed in terms of 
the effective-range parameters by matching to the effective-range expansion for low-energy 
fermion-fermion scattering ^36']: 



Co = ^, C, = Co^, and C^ = i^, (44) 

where (op) is the s-wave (p-wave) scattering length and is the s-wave effective range, 
respectively. If the effective range parameters are all of order the interaction range, then 
the EFT is said to be natural. In a uniform system at finite density, the mean inter- 
particle spacing tq provides a length scale for comparison; the ratio a^/ro ~ kptts provides a 
dimensionless measure of density, where kp is the Fermi momentum. In the dilute regime, 
kpCis 1 serves as an expansion parameter, as realized by the EFT of Ref. where the 
diagrams of Fig. QJeach contributed to precisely one order in the energy density: LO or Si 
is C»(4), NLO or £2 is 0{kl), and NNLO or £3 is O(fcf). 

In extending these results to a finite system, the simplest approximation that can be 
invoked is the local density approximation (LDA). In Ref. |33|, the results for a uniform 
system were used in the LDA to evaluate the F/ and hence the energy through NNLO. We 
simply quote the results here. The LO diagram [Fig. ^a)] is the Hartree-Fock contribution, 
which is purely local. Thus, it is an exact evaluation and Eq. (|34j) can be used directly for 
the energy contribution: 

Ehf[p{^)] = I d'^ £M^)] = Co I rf^x |p(x)p . (45) 

The contributions to the energy from NLO and NNLO diagrams are computed in LDA by 
simply integrating the corresponding uniform energy densities evaluated at the local density 
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0: 

+ (62 r, + hai + h al) ^ J rf^x [p{^)f' , (46) 
where the dimensionless bk are 

IOtt \ ^ / 

(e 2 \ 5/3 
-^j (^0.0755 (z/ - 1) + 0.0574 (z/- l)(z/- 3)) . (47) 

The numerical constants in the last line of Eq. (jTfjl were obtained by Monte Carlo integration 
0. 



C. Including r in Hartree-Fock Diagrams 

The contributions to i?^°^[p(x)] from the Hartree-Fock graphs containing the C2 and C2 
vertices [Figs, ^i) and (j)], which have gradients, were approximately evaluated in the LDA 
in Eq. (j46j) to obtain a functional of p alone. In the DFT formalism generalized to include 
the kinetic energy density t(x), however, we observe that they can be evaluated exactly in 
terms of p(x) and r(x) (for closed shells). These contributions to M/3[Jo)''?o] can be simply 
expressed in terms of gradients acting on the Kohn-Sham Green's functions [cf. Eq. fl32|) ]: 

C2 z/ (z/ - 1) {(Vi - V2)' + (V3 - V4)'} 
+ 2C!, z/ (z/ + 1) {(Vi - V2) • (V3 - V4)} 

+ ■ ■ ■ (48) 

Equation (PH|) is evaluated by carrying out the gradients and then setting all of the Xi equal 
to X. For this purpose, the replacement 

occ. 

z^^fc(x)^*(x') (49) 

k 



W3[Jo,Vo] = -T^ 



(i^X 



^^(^4, 4)^^(^3,0;+) 



Glix,x') 
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can be made and, since in the end we integrate over x, we can partially integrate any term. 
It is not difficult to write the resulting expressions in terms of p(x), r(x), and a current 
density 



occ. 



j(x) = -z/5^[V^^(x)^,(x) -^^(x)VV^,(x)] . (50) 



2 

k 



The current density vanishes when summed over closed shells, leaving simple expressions for 
the Hartree-Fock energy functionals: 

i?cjp(x), r(x)] = j rf^x [p(x)r(x) + \ {V pf] , (51) 



[p(x), r(x)] = ^jd'^ [p(x)r(x) - \ {V pf] . (52) 
The dimensionless constants B2 and i?3 are given by 

B, = ,(^, B3 = 2.i^. (53) 

The total contribution to the energy from NLO and NNLO diagrams is 

E,[p(x),r(x)] = H |rf=^x[p(x)]^/=^ + i?cJp(x),r(x)] 

+ E^.[p(x),r(x)] + ^/ c^^x [p(x)]«/3 . (54) 

Since this functional now has the semi-local r(x) as one of its ingredients, it represents a step 
beyond LDA, even though the contributions from Figs, mb), (g), and (h) are still evaluated 
with the LDA prescription. The next step in the DFT/EFT program will be to develop 
systematic expansions to these contributions in terms of p and r. 

For now, we carry out the DFT/EFT formalism in the effective action framework using 
the hybrid functional with r(x). The full effective action is given by 

00 

f [p, r] = f o[p, r] + f i[p, r] + ^ f ,[p, r] . (55) 

fe=2 

We proceed to calculate the sources using Eqs. (j39j) . (pUj) . and (j54j) to NNLO; ffist Jo(x) is 

Jo(x) = ^(^fi[p,r] + |^f,[p,r]j =-^(^HF[p]+i?c[p,r]) 

= 1^ p(x) - I 6, ^ [p(x)]^/3 - I ^[p(x)]V3 

- (i?2 r, + 53 aj) ^ r(x) + (SE^ r, - B, aj) ^ VV(x) , (56) 

and then 77o(x) is: 

"^'^""^ = ^^"^['"^ + ^'^f'"' ""^^ = " + ""^^ 2^ P^""^ • ^^'^^ 
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The spatially dependent effective mass M*(x) is therefore 



1 



2M* X 



2M 



2M 



4z/ 



Co 



4i/ 



P(x) 



(58) 



where Eq. ()44|) has been used. An expression for the total binding energy (through NNLO) 
follows by substituting for Jo(x) and r7o(x) in Eq. and then using Eq. and Eq. (fT3j). 



E[p(x),r(x)] 



occ. 



\ iu — \) An as r , ,-,9 4 , 



[P(x)] 



7/3 



+ ^4 ^ [/^(^)]'^' + (3^2 - aj) ^ [Vp(x)] 



8M 



+ (52a^r, + 53a3)^p(x)r(x)} 



(59) 



The noninteracting case (Cj = 0) uses Jo(x) = 0, 77o(x) = and the first term in Eq. (j^HI)- 
Leading order (LO) uses the first term in Eq. (|56|) and the first two terms in Eq. (j59|) . and 
so on for NLO and NNLO. An alternative expression for the energy is obtained by using the 
second part of Eq. ^ followed by Eqs. (IHIH) . and (ITK|) : 



E[p(x),r(x)] 



— r(x)+.(x)p(x) + -L_J__[p(x)P 



(i?2 al r, + ^3 aj) ^ p(x) r(x) + (3i?2 r, - i?3 4) ^ [Vp(x)]2 



8M 



8/3 1 



(60) 



D. Comparison to Skyrme Hartree-Fock 

A nucleus is a self-bound system, so the external potential f (x) =0. To compare the 
DFT/EFT functional to the conventional Skyrme energy density functional of Eq. (jSJ, we 
set the spin multiplicity v = A and use Eq. to rewrite Eq. (jUUI) in terms of the Cj's, 
obtaining the energy density: 

£:[p(x), r(x)] = ^r(x) + \c, [p(x)]2 + ^(SC^ + 5C;)p(x) r(x) + ^(gC^ - 5C;)(Vp)2 

+ ^l||[/^W^/^ + &4^[p(x)]«/^ (61) 

We observe that we get all the terms of £^5/^(x) from Eq. (0) except for the one with coefficient 
^3 and the spin-dependent terms, if we make the correspondence to ^ C'o; ^ and 
^2 ^ C^- This correspondence is not surprising since the Skyrme interaction was originally 
motivated as a low- momentum expansion of the G matrix j2^. The two additional terms 
in Eq. (j?)T|) of the form p(x)^^" come from correlations (i.e., terms beyond Hartree-Fock), 
but there is no direct association with the ^3 term in the Skyrme energy density, which 
was originally motivated as a three-body contribution (so a = 1). However, it is clear that 
the Skyrme functional is incomplete as an expansion; a direct connection to microscopic 
interactions by matching to an EFT will include at least these additional terms. 



14 



The generahzation of the DFT/EFT to include spin and isospin dependence is straight- 
forward. If one writes a complete set of four-fermion terms with cr and r matrices in the 
EFT Lagrangian, there are redundant terms because of Fermi statistics. In conventional 
discussions of the Skyrme approach, this observation is typically cast in terms of antisym- 
metrization of the interaction jisl . [l9l |. For a path integral formulation of the DFT/EFT, 
using Fierz rearrangement is a convenient alternative. We illustrate the procedure for the 
leading spin dependence. 

First consider just spin-1/2 (no isospin, so v = 2). We expand the product of Grassmann 
fields V'jV^j (^ aiid j are spin indices) in the complete basis of 5ij and af^- (a = {1,2,3}), 
identifying the coefficients by contracting in turn with 5ij and a^j. The result, with minus 
signs from interchanging Grassmann fields, is 

M = - • (62) 

If we substitute this result into {ip'^ipY = ipl{tpitpj}tpj, we find 

(V^W = -^(V^W-^(^W)% (63) 

or 

(We could also start with {ipja1i,4'i'){'ipjcrjji'ipj') and obtain the same result with a bit more 
effort). Therefore 

Co(^V)' + C^ii^^fTijf = (Co - 3Co^)(V^V)2 (65) 

and the single term Co{ip'^ipy with Cq = Cq — 3Cq yields the same results for all diagrams 
as the original two terms. 

For the z/ = 4 case with spin and isospin, we perform a similar procedure to find 

^^c.i^]p = -\{^P^mAp - \{^^cr^^K6^p - ^(^VV)5.,r^^ - i(^VvV)a,^r^;3 • (66) 

Substituting into any two of {ip'^ipy, {ip'^cnpY, {ipWipY, and {ip'^CTTipy, we find two inde- 
pendent relations, which can be solved simultaneously to find: 

{^^T^f = -2(V^V)^ - (^V^)2 (67) 
{^^(TT^'f = -3(V^V)^ , (68) 

which allow us to eliminate explicit dependence on r matrices in favor of just two indepen- 
dent couplings: 

Co(v^V)' + c^{^^^^? + Co^(^v^)2 + c^^i^^cTT^y 

= {Co - 2Cl - 3Co^-)(V^V)' + (Co" - CDi^lj^cT^jf . (69) 

This agre es ( of course) with the usual discussion in terms of antisymmetrized interactions 
(e.g., see or ^^). The choice to eliminate the r terms is purely conventional. The 
convention with Skyrme interactions is to choose the independent couplings to be ti and Xi, 
which multiply terms in the effective interaction in the combination ti(l + XiPa), with 
the spin exchange operator. Extending to more derivatives and the spin-orbit (and tensor) 
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terms follows systematically in the EFT approach, but introduces many more constants; 
complete sets of contact terms with more derivatives in the EFT Lagrangian can be found in 
Refs. (s^l and js^. The prohferation of constants leads to a clash of philosophies between the 
minimalist, phenomenological approach (use as few terms as possible), which is necessarily 
model dependent, and the model-independent EFT approach (use a complete set of terms). 

The similarities of the successful Skyrme functional and the DFT functional for a dilute 
Fermi gas prompts an analysis of typical Skyrme parameters as effective range parameters. 
One can use the association of the tj's and the CiS along with numerical values from suc- 
cessful Skyrme parameterizations (e.g., Ref. j2l[) to estimate "equivalent" values of a^, r^, 
and ttp. We find that ^ — 2-3fm, which is about the inverse pion mass and is much 
smaller than the large, fine-tuned values of the free-space nucleon-nucleon interaction (how- 
ever, some Skyrme parameterization have an "equivalent" of 5fm or larger). However, 
kptts is still significantly larger than unity inside a nucleus, which precludes a perturbative 
dilute expansion. Interestingly, the equivalent and values have magnitudes consistent 
with what one might expect from the nuclear hard-core radius (with ap < 0), leading to kpTg 
and kpap less than unity. 

III. RESULTS FOR DILUTE FERMI SYSTEM IN A TRAP 

In this section, we present numerical results for the dilute Fermi system comprised of a 
small number of fermions confined in a harmonic oscillator trap. We compare the nature of 
the convergence of the EFT in a finite system both qualitatively and quantitatively to the 
analysis done purely in the LDA js^, which means the effect of treating the Hartree-Fock 
contributions at NNLO [Fig. HJ^i) and (j)] exactly. 

A. Kohn-Sham Self-Consistent Procedure 

We restrict our calculations to finding the Kohn-Sham orbitals for closed shells, so the 
density and potentials are functions only of the radial coordinate r = |x|. Note that the 
basic procedure is the same one used for closed-shell nuclei in Skyrme-Hartree-Fock jisj ]. 
even though the DFT/EFT can include correlations to any order. The Kohn-Sham iteration 
procedure is as follows: 

1. Start by solving the Schrodinger equation with the external potential profile v{r) for 
the lowest A states (including degeneracies) to find a set of orbitals and Kohn-Sham 
eigenvalues {ipky^k}- 

2. Compute the density and kinetic energy density from the orbitals: 



A 




(70) 



k=l 




k=l 



(71) 
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3. Using Eqs. find Jo{r) and rio{r). Evaluate the local single-particle potential 

Vs[p{r), r(r)] = Vs{r) = v{r) - Jo{r) (72) 
at the chosen level of approximation (e.g., NLO) and the "effective" mass : 

4. Solve the Skyrme-type Schrodinger equation for the lowest A states (including degen- 
eracies), to find {ipky^k} as before: 

[-V^^^^V + Vs{r)] M^) = ekM^) ■ (74) 

5. Repeat steps 2.-4. until changes are acceptably small ("self-consistency"). In practice, 
the changes in the density are "damped" by using a weighted average of the densities 
from the {n — l)th and nth iterations: 

p(r)=/3p„_i(r) + (l-/5)p„(r) , (75) 

with 0< f3 <1. 

This procedure has been implemented for dilute fermions in a trap using two different 
methods for carrying out step 4. The Kohn-Sham single-particle equations are solved in 
one approach by direct integration of the differential equations and in the other approach 
by diagonalization of the single-particle Hamiltonian in a truncated basis of unperturbed 
harmonic oscillator wavef unctions. The same results are obtained to high accuracy. [Note 
that other methods used for Skyrme-type equations, such as the conjugate-gradient method 
(s^ . can also be directly applied.] 



B. Fermions in a Harmonic Trap 

The interaction through NNLO is specified in terms of the three effective range parameters 
ttg, Tg, and ttp. For the numerical calculations presented here, we consider the natural case 
of hard-sphere repulsion with radius R, in which case ag = ap = R and = 2R/3, and also 
the case with = 2as, so we can emphasize the effect of M*{r)/M significantly less than 
unity. 

Lengths are measured in units of the oscillator parameter b = a/ h/Muo, masses in terms 
of the fermion mass M, and h = 1. In these units, hut for the oscillator is unity and the 
Fermi energy of a non-interacting gas with filled shells up to Np is Ep = [Np + 3/2). The 
total number of fermions A is related to Np by 

A = ^{Np + l){Np + 2){Np + 3) . (76) 

Since we have only considered spin-independent interactions, our results are independent of 
whether the spin degeneracy u actually originates from spin, isospin, or some fiavor index. 

With interactions included, single-particle states are labeled by a radial quantum number 
n, an orbital angular momentum / with z-component rui, and the spin projection. The radial 
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TABLE I: Energies per particle, averages of the local Fermi momentum kp, and rms radii for 
sample parameters and particle numbers for a dilute Fermi gas in a harmonic trap. See the text 
for a description of units. The scattering length is fixed at Ug = 0.16 and the effective range is set 
to rg = 2as/2> when Op ^ 0. Results with the DFT functional including r are marked "r-NNLO." 



V 


Nf 


A 


Up 


E/A 


{kp) 


7F) 


approximation 


2 


7 


240 




7.36 


3.08 


2.76 


LO 


2 


7 


240 




7.51 


3.03 


2.81 


NLO (LDA) 


2 


7 


240 


0.00 


7.52 


3.02 


2.82 


NNLO (LDA) 


2 


7 


240 


0.16 


7.66 


2.97 


2.87 


NNLO (LDA) 


2 


7 


240 


0.16 


7.65 


2.97 


2.87 


r-NNLO (LDA) 


2 


7 


240 


0.32 


8.33 


2.76 


3.10 


NNLO (LDA) 


2 


7 


240 


0.32 


8.30 


2.77 


3.09 


T-NNLO (LDA) 



functions depend only on n and /, so the degeneracy of each level is z/x (2/ + 1). The solutions 
take the form (times a spinor, which is suppressed) 



Unl(r 



where the radial function Uni{r) satisfies 
1 (P dr]Q f 1 d 



2M*{r)dr'^ dr \r dr 
and the Uni^s are normalized according to 



+ v.Jr) + 



/(/ + !) 
2M*(r)r2 



\uniir)\'^dr = 1 . 



Thus the density is given by 



nl 



and the kinetic energy density is given by |5; 



nl 



dR, 



nl 



dr 



+ 



1(1 + 1) 



Rniir)\' 



(77) 



(78) 



(79) 



(80) 



(81) 



The interactions are sufficiently weak that the occupied states are in one-to-one correspon- 
dence with those occupied in the non-interacting harmonic oscillator potential. 



C. Numerical Results 

Here we compare density distributions at zero temperature for the LDA analysis [3c 
to those from evaluating Figs, ^i) and (j) exactly at NNLO (where the differences first 
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FIG. 2: NNLO Kohn-Sham density distributions for a dilute gas of fermions in a harmonic trap 
with degeneracy u = 2 fihed up to N'p = 7, which imphes there are 240 particles in the trap. The 
scattering length is = 0.16 and the effective range is = 2as/3. Results for two values of Up 
are compared for the LDA /9-only functional (p-DFT) and with r (pr-DFT). 
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FIG. 3: Deviation of pr-DFT from p-DFT results at NNLO for the same systems as in Fig. |^ 
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FIG. 4: NNLO Kolin-Sham kinetic-energy-density distributions for the same systems as in Fig. |21 
The upper three curves are for ap = as and the lower three curves are for Op = 2as. 



appear in the present analysis). This comparison obviously makes sense only if at least 
one of Ts and ap is non-zero. We choose 240 particles as a representative example (other 
numbers of particles give qualitatively similar results). Densities at different orders in the 
DFT expansion were shown in Ref. j33|. Results for the energy per particle E/A, average 
Fermi momentum (k-p) and the rms radius are given in Table HI 

In Fig. 121 we compare the densities at NNLO in the Kohn-Sham formalism with p-only 
functionals (the LDA calculation from Ref. js^]) to the same system with the pr functionals. 
For hard-sphere scattering (for which ap = as), the density curves are almost indistinguish- 
able. If we plot the difference on an expanded scale (see Fig. El), we can see a small amplitude 
oscillation. The close agreement is not surprising given that the source of the difference is 
the NNLO Hartree-Fock terms, so the difference itself is higher order in the EFT expansion 
(note that the p-DFT and pr-DFT NNLO energies in Table jll differ only by 0.01). We can 
magnify the difference by considering ap = 2as, which multiplies the corresponding Hartree- 
Fock term by a factor of eight (which also implies that the coefficient is unnaturally 
large). For this case, the difference in Fig. El is visible and significant oscillations are seen in 
Fig. 121 The increased oscillation is analogous to the difference between Thomas-Fermi and 
Kohn-Sham DFT densities (shown in Ref. [33]), although not as dramatic. The explanation 
is also analogous: the pr-DFT captures more non-locality into the DFT functional. 

The kinetic energy densities for these cases are shown in Fig. 01 When calculated from 
the Kohn-Sham wave functions, r(r) is quite similar for the p-only and pr calculations. Also 
shown in this figure is the leading contribution from the semiclassical approximation, which 
reproduces the "exact" kinetic energy densities except near the origin. 

Only in the pr case is the effective mass M*(r) different from unity; it is shown in Fig. jJl 
for the two values of ap. The values close to the origin are in the range of those obtained in 
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FIG. 5: NNLO Kohn-Sham effective mass distributions for the same systems as in Fig. |2I Results 



for two values of ap are shown. 



Skyrme functionals fit to nuclear data. In those functionals, the value of M* is associated 
with the single-particle energy levels. This correspondence is seen in the single-particle 
energy spectra in Fig. IHl If we take the uniform limit (with the external potential turned 
off), the single-particle energies for momentum k differ according to 



2M 



(82) 



Thus, for positive and ap, the pr levels will always lie lower except at the Fermi surface 
(where they must be equal). 

This comparison demonstrates how the Kohn-Sham formalism can be misinterpreted. 
Even though the single-particle energies differ significantly, they are not observables. Indeed, 
the true bulk observables calculated in the DFT framework, those in Table IH are barely 
distinguishable. One could ask whether the single-particle levels in some representation are 
"better" than in other representations. In particular, they can be compared to the energy 
spectrum corresponding to the poles of the exact Green's function, which can be constructed 



in terms of the Kohn-Sham Green's function |48|, jH^. If this construction is carried out in 
the present approximations for both the p-only and the pr formalisms, the pr single-particle 
spectrum is obtained in both cases js^. It is not clear from the present calculation how 
close the full and Kohn-Sham spectra would be if the LDA were relaxed [e.g., for Fig. ^b)]. 
It would be useful to add spin-orbit interactions and then to study Kohn-Sham spin-orbit 
splittings near the Fermi surface, since such splittings are commonly fit in mean-field models. 
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FIG. 6: Comparison of selected single-particle energy spectra for the same systems as in Fig. |2j 



D. Power Counting and Convergence 

The effective field theory approach allows us to estimate contributions to the energy. 
At each successive order in the EFT expansion, the low-energy constants (LEC's) can be 
estimated using naive dimensional analysis, or NDA [s^]. In the case of a short-range force 
with a natural scattering length, the underlying momentum scale A ~ where R is the 
range of the potential, is the basic ingredient in the NDA. The estimate of two-body Hartree- 
Fock energy contributions from a given term in the Lagrangian can be found by replacing 
ip'^ip by the average density (and including an appropriate spin factor) and the coefficient 
by the natural estimate ~ 47r/MA^*+^ [3^. As an example, the NDA estimate of the 
Hartree-Fock energy per particle at LO was computed from Eq. (j45j) as : 



E 



HF 



A 



NDA 



2 V 



An 
MA 



(P(x)) , 



(83) 



with A = 1/ R = l/ag for a hard-sphere potential. The Thomas- Fermi result can be used 
to find the average density, or one might just take the actual computed average value (the 
results will differ by much less than the uncertainty in the estimate). The NLO estimate was 
found by multiplying the LO result by {k-p)as, where (kp) is the average Fermi momentum 
j^^l- At NNLO, we have three terms. The p^^^ LDA term was estimated by multiplying the 
LO result by((A;F)as)^, and the other two terms (pr and Vp) arising from Hartree-Fock at 
that order was estimated directly as in the case of LO (using (pr) ^ (p)(r)). 

In Figs. [7| and |H1 estimates and actual contributions are shown for u = 4, A = 140, 
tts = 0.10 (for which {kpa^) ^ 0.24) and for z/ = 2, A = 240, = 0.16 (for which {kpa^) ^ 
0.5). Square symbols denote estimates based on naive dimensional analysis, with error 
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bars indicating a 1/2 to 2 uncertainty in the estimate. Actual contributions to the energy 
per particle from each of the orders are shown as round symbols, i.e., the actual NLO 
contribution is |-Enlo — -^loIM- At NNLO, we plot estimated contributions from the p^^^ 
LDA, p r, and gradient terms separately. The latter gives a very small contribution consistent 
with its NDA estimate, and it has been multiplied by ten in Figs. |7| and |H1 to fit them on 
the graphs. The two sets of estimates and results correspond to Op = (hard sphere) and 
Op = 0. 

We see from Figs. |71 and |H1 that the actual results for LO, NLO and NNLO agree well 
with NDA estimates, including the new gradient contributions, with one exception. The 
exception is the p^^^ LDA estimate for the z/ = 2 system, which greatly overestimates the 
actual contribution at that order due to an accidental cancellation when v = 2 between 
the two terms in the 64 coefficient in Eq. (jTTj) [33]. In general, however, we can use these 
estimates to reliably predict the uncertainty in the energy per particle from higher orders. 

While it's clear that nuclei are not perturbative dilute Fermi systems with natural free- 
space scattering lengths, there is phenomenological evidence that power counting can apply 
to energy functionals that are fit to bulk nuclear properties. Phenomenologically successful 
functionals, both of the Skyrme type and covariant, have Hartree terms that are consistent 
with Georgi-Manohar NDA for a chiral low-energy theory 5^ 57, 5^. This involves power 
counting with two scales: the pion decay constant /^r and an underlying scale for short- 
range physics A, which empirically (for these functionals) is around 600 MeV. (Equivalently, 
there is an additional large dimensionless coupling g A/ f^^ that enters in a well-prescribed 
manner). Thus the functionals do take the form of a density expansion (with parameter 
p/ f^A) as well as a gradient expansion, with the same hierarchy as illustrated here. A major 
goal of future investigations will be to elucidate the nature of the density expansion for finite 
nuclei and to connect it to the underlying chiral EFT. 



IV. SUMMARY 

In this paper, the EFT-based Kohn-Sham density functional for a confined, dilute Fermi 
gas was extended by incorporating the kinetic-energy-density r(x) into the formalism. The 
generating functional is constructed by including, in addition to a source J(x) coupled to 
the composite density operator ip'^ip, another source //(x) coupled to the (semi-local) kinetic 
energy density operator Vip"^ ■ Vip. A functional Legendre transformation with respect to 
the sources yields an effective action of the kinetic energy density r as well as the fermion 
density p. This construction also serves as a prototype for including additional densities 
currents (such as separate proton and neutron densities or a spin-orbit current). 

This extension sets the stage for the construction of energy functionals and Kohn-Sham 
equations that go systematically beyond the local density approximation (LDA). As a first 
step, we included the exact Hartree-Fock (HF) contribution at NNLO for a natural dilute 
Fermi system but treated non-HF contributions in LDA. This exact HF contribution provides 
explicit dependence on r and on the gradient of the density, unlike the LDA. The EFT 
expansion for a confined, finite system with natural effective range parameters is controlled 
by those parameters (a^, r^, ap, . . . ) times an average Fermi momentum and by the gradient 
of the density. Thus the expansion is perturbative in the sense of being a density expansion 
but is not perturbative in an underlying potential. (This is clear since our prototype system 
is hard spheres, which yields infinities even in first-order perturbation theory!) An error 
plot of contributions to the energy per particle versus the order of the calculation showed 
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FIG. 7: Energy estimates (squares and error bars) for = 4, Cs = 0.10, A = 140 particles compared 
with actual values (circles) for a hard-sphere gas (sohd) and with the p-wave scattering length equal 
to zero (shaded). 




T-NNLO 

FIG. 8: Energy estimates (squares and error bars) for u = 2, = 0.16, A = 240 particles compared 
with actual values (circles) for a hard-sphere gas (solid) and with the p-wave scattering length equal 
to zero (shaded). 
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that we can reliably estimate the truncation error in a finite system, including the gradient 
terms. 

The ground-state energy functional and Kohn-Sham single-particle equations constructed 
here take the same form as those in Skyrme Hartree-Fock calculations (not including the 
spin-orbit contribution, which can be added with a similar generalization). There is work 
supporting the Skyrme (and also covariant 5^|) mean- field energy functionals as density 
expansions. Low-energy effective theories of QCD are expected to have a type of power 
counting. As shown in Ref. js^, power counting in the energy functional based on chiral 
naive dimensional analysis is fully consistent with phenomenological parametrizations. The 
value of including the kinetic energy density as a functional variable is not conclusive in the 
present example, but it does lead to a single-particle spectrum closer to that of the exact 
Green's function. In general, summing more nonlocalities without decreasing the numerical 
efficiency of the calculation is a plus. 

Even in the simplest low-density expansion, there are contributions at all powers of the 
Fermi momentum, which means fractional powers of the density. In the Skyrme parametriza- 
tion, there is only one such term. This is presumably a balance between phenomenological 
accuracy and a desire to maximize predictive power through limiting parameters. For repro- 
ducing bulk properties of stable nuclei, the Skyrme functional is only sensitive to a relatively 
small window in density, which allows significant freedom for a phenomenological parameter- 
ization. The future challenge will be to see if additional terms motivated by power counting, 
which may become important for extrapolating far from stability, can be determined. 

The immediate next steps are to develop a derivative expansion to evaluate diagrams 
such as the NLO "beachball" diagram beyond LDA and to test it for convergence, and to 
generalize the DFT formalism to include pairing. Work is in progress on these extensions, 
which will be directly relevant for nuclear applications j^. In addition, in order to adapt 
the density functional procedure to chiral effective field theories with explicit pions, we will 
need to extend the discussion to include long-range forces. Kaiser, Weise, and collaborators 
have already generated functionals of the Skyrme type using free-space chiral perturbation 
theory together with the density matrix expansion of Negele et al. jH^. It is not clear that 
the power counting is consistent in those calculations, but future development along these 
lines is certainly warranted. Finally, a systematic solution to the large scattering length 
problem for trapped atoms also remains a challenge. 
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